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ABSTRACT 

Context. Magnetic neutral points are potential locations for energy conversion in the solar corona. 2D X-points have been widely 
studied in the past, but only a few of those studies have taken finite plasma beta effects into consideration, and none of them look at 
the dynamical evolution of the system. At the moment there exists no description of the formation of a non-force-free equilibrium 
around a two-dimensional X-point. 

Aims. Our aim is to provide a valid magnetohydrostatic equilibrium from the collapse of a 2D X-point in the presence of a finite 
plasma pressure, in which the current density is not simply concentrated in an infinitesimally thin, one-dimensional current sheet, as 
found in force-free solutions. In particular, we wish to determine if a finite pressure current sheet will still involve a singular current, 
and if so, what is the nature of the singularity. 

Methods. We use a full MHD code, with the resistivity set to zero, so that reconnection is not allowed, to run a series of experiments 
in which an X-point is perturbed and then is allowed to relax towards an equilibrium, via real, viscous damping forces. Changes to 
the magnitude of the perturbation and the initial plasma pressure are investigated systematically. 

Results. The final state found in our experiments is a "quasi-static" equilibrium where the viscous relaxation has completely ended, 
but the peak current density at the null increases very slowly following an asymptotic regime towards an infinite time singularity. 
Using a high grid resolution allows us to resolve the current structures in this state both in width and length. In comparison with the 
well known pressureless studies, the system does not evolve towards a thin current sheet, but concentrates the current at the null and 
the separatrices. The growth rate of the singularity is found to be f", with < D < 1. This rate depends directly on the initial plasma 
pressure, and decreases as the pressure is increased. At the end of our study, we present an analytical description of the system in 
a quasi-static non-singular equilibrium at a given time, in which a finite thick current layer has formed at the null. The dynamical 
evolution of the system and the dependence of the final state on the initial plasma and magnetic quantities is discussed, as are the 
energetic consequences. 
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1. Introduction 

Understanding the nature of hydromagnetic equilibria is a key 
issue for the study of solar magnetic fields such as those of the 
solar corona. Magnetic relaxation and field extrapolation are top- 
ics which clearly involve magnetic equilibria, and many other 
areas of study start from and/or aim to evolve to an equilibrium 
magnetic field, namely, magnetoacoustic wave propagation and 
dissipation, magnetic reconnection, or current sheet formation. 
However, the majority of these studies simply involve force-free 
equilibria, i.e. equilibria in which only magnetic forces are con- 
sidered and forces due to plasma pressure, gravity, or velocity 
are all ignored. 

In general, such assumptions are probably very reasonable 
in the majority of the low beta solar corona. However, in certain 
regions, such as near null points (where |B| is very small) and 
in the chromosphere, the effects of the plasma pressure become 
important. A study of the non-force-free relaxation of unbraided 
magnetic fields has started to be studi ed numerically and ana- 
lytica lly in the first paper of this series dFuentes-Fernandez et al.l 

l2oToh . 

Magnetic energy release in the solar corona has for long 
been, and still remains, a key unresolved problem in solar 
physics. Magnetic null points play an important role in this re- 
lease of magnetic energy, since they are likely to evolve by col- 
lapsing towards a tangential field discontinuity with a strong 



current, eventually leading to dissipation via magnetic recon- 
nection. Such processes have been widely studied both analyt- 
ically and numerically, with direct applicatio ns to solar envi- 
ronments such as in the CME breakout model (lAntiochos et al.l 
1999), which has been applied extensively in the last decade 
(e.g. For bes et al.ll2006t IZuccarello et al.ll2009l) . and in other in- 
terplanetary scen arios such as the reconnec tion site in the Earth's 
magnetotail (e.g. lHesse & Schindled200"Tl) . Also, they have been 
used in wave propagation experime nts involving a zero beta 
plasma McLaughlin & Hood (2004) and a finite beta plasma 
( McLaughlin & Hood 2006) . finding in both cases that the waves 
wrap around the null point, causing an exponential build up of 
current density at the location of the null. On the other hand, 
2D rec onnection has been studied for decades st arting with 
Dungev (195 3[) and followed up by rnany (e .g. Parker 19571 
Sweet, 1958.: IPetschekl [1961 iBiskamd Il986t jPriest & ForbesI 



il986t]craig & Hentonlll994HShibata et al.lll993i etc.) 



Furthermore, null points have been found to have a rea- 
sonably high population density in the solar corona, by 
iLongcope & Parnelil(l2009l) . In two dimensions, in particular, the 
topology of the current sheet formed after an X-type null point 
collapse has been well studied analytically for potential and 
force- free fields by several autho rs (e.g. Dun gev 1953: Greera 
1 965t ISomov & Svrovatskiill 19761; ICraidll994l;lBungev & Priesfl 
.199.5h . 



2 



Jorge Fuentes-Fernandez, Clare E. Pamell and Alan W. Hood: MHD Dynamical Relaxation of Coronal Magnetic Fields 



The aim of this paper is to provide a valid magnetohydro- 
static equilibrium from the collapse of a two-dimensional X- 
point, which are simpler scenarios than 3D null points, although 
the latter will be studied in a follow-up paper. Under ideal, non- 
resistive conditions, the energy bound up in the global magnetic 
field has to manifest itself as localized accumulations of current 
density. 

It is well known that under the cold plasma approximation 
(e.g. zero plasma beta), an initially perturbed X-point field re- 
laxes to a potential equilibrium with a Y-type infinitesimally thin 
current sheet where the current is zero everywhere except within 
the magnetic tangential discontinuity, where it develops a singu- 
larity of the form /, - S{A. -A^q). These potential configuration s 
are de scribed b y Green ( 1965| and lSomov & Svrovatskiil(ll976h . 
Later. iBungev & Priesti (|199 5) expanded these solutions for po- 
tential and force-free fields giving a general expression for these 
force-free current sheets. Latter studies have found the forma- 
tion of localised i nfinite current layers in the Earth's magnetotail 
dBirn et al.ll2003h . relevant for the initiation of the subsequent 
energy release phase. 

fpriedel et al. ( 1996) and iGrauer & MarUanil (Il998h sfiidied 
two-dimensional current singularities in incompressible flows 
using an Adaptive Mesh Refinement (AMR) code which per- 
mitted them to obtain very highly localised resolutions at the lo- 
cations of maximum current. They found an exponential growth 
for the current at the null p oint. 

An analytical result by ( lKlapperll997h has rigorously shown 
that, in 2D ideal incompressible plasmas, a singularity of the 
current density will take an infinite amount of time to develop, 
unless driven by a pressure singularity occurring outside a neigh- 
borhood of the null point. In this paper we aim to show numeri- 
cally that this is also true for compressible plasmas with a finite 
plasma pressure. 

First evidence of current sheets extending along the sep- 
aratrices in shear e d mag netic field structures were studied by 
IZwingmann et al.l (1 1985 1) for force-free equilibria, where they 
found mathematical singularities in the current sheet which they 
interpreted as terms that "would become large in a real physical 
situation". Later, Vekstein & Priest ( 1993) made a mathematical 
analysis of the magnetic field around cusp-points, after the shear- 
ing of a magnetic field with an X-type null point, and suggested 
a form of a negative power law for the resulting singular current 
de nsity as a fu nction of A^. 

lRastatter~ et al. ( 1994|) considered for the first time the efifects 
of pressure perturbations in numerical experiments on the ideal 
relaxation of two-dimensional magnetic X-points, and studied 
the development of current layers with singular current densities 
in which, in the relaxed state, the initial X-point was replaced 
by either a T-point or a cusp-point geometry. For the relaxation, 
they used a frictional code, which damped the kinetic energy out 
of the system by adding a fictitious relaxing term to the momen- 
tum equation of the form -a:v, but without any associated heating 
term in the energy equation. Their X-point relaxed to a singular 
equilibrium with a plasma pressure jump across the separatrices, 
and with current layers extending along the separatrices. They 
argued that the finite width of their current sheet was due to the 
finite difference method in their numerical approach rather than 
being real, but found, nevertheless, the integrated current den- 
sity over the sheet width (referred as to surface current) to be 
constant a long each whole separatrix . 

Later, 'Craig & Litv inenkol (l2005h reconsidered the problem 
of the relaxation of two-dimensional magnetic X-points and the 
formation of current singularities in non-force-free equilibria, 
and evaluated the strength of the current singularity at the end 



of their relaxation. Again, they made use of a frictional code 
with a fictitious damping term, -/cv, added to the momentum 
equation, but with no heating term in the energy equation, as- 
suming the polytropic model p ~ p^, which imposes a condi- 
tio n of adiabaticity to th e process. In analogy with the results 
of iRastatter et al.l (Il994t) . they found a distribution of current 
density extended along the magnetic separatrices, which they 
claimed to be almost uniform. They repeated their study for a 
number of different grid resolutions, and presented a logarith- 
mic increase of the peak current with the number of grid points, 
at the same time as the area of the current layer above a given 
value for the current showed a logarithmic decrease. Hence, the 
current layer itself became narrower with higher resolution. 

Then, they evaluated the scalings of the peak current for dif- 
ferent values of the background plasma pressure of the system, 
finding a weakening of the growth of the peak current density as 
the plasma pressure was enhanced. That is, a singularity is harder 
to achieve the higher the value of the plasma pressure, although 
the presence of a non-zero plas ma beta would not prevent a sin- 
gularity forming. The results of^Craig & Li tvinenkol (l2005h were 
fully reproduced by Pontin & Craig (2005), as part of their study 
of current singularities at 2D and 3D nulls, using a Lagrangian 
code, instead of the Eulerian code used in ICraig & Litvinenkol 
(I2OOI . 

Here, we study the dynamical evolution of magnetic X- 
points under the frozen-in assumption, which is driven by vis- 
cous forces which involve a heating term that ensures energy 
conservation. The forces of gravity are neglected, as the pres- 
sure scale-heights in the solar corona (of the order of lOOMm) 
are far larger than our length-scales, here understood (given our 
line-tied boundaries) as the height of the possible magnetic null 
points above the base of the corona, whic h might be of the order 
of IMm (see iLongcope & Parnelll 120091) . We show how a type 
of singularity wants to be formed in the non-force-free case, in 
agreement with the numerical studies of Rastatter et al.l (IT99I 
and Craig & Litv inenko ( 2005). Our numerical results show how 
the initial X-point collapses to a cusp-like geometry in which the 
current density accumulates around the neutral point and along 
the four separatrices. Again, the results agree, in this aspect, 
with the previous numerical works of non-force-free X-point 
collapse. However, we attempt to go a step further by looking at 
the time evolution of the field, running a series of very high res- 
olution experiments, which allow us to look closer at the current 
accumulations. Also, following the lines of Vekstein & Prie^ 
(II993I) . we propose a mathematical form for the final equilib- 
rium state, involving a finite current layer ending in cusp-points 
at both extremes. 

The paper is structured as follows: In Sec.|2| we present the 
initial setup and the details of the numerical experiments. In Sec. 
|3]we state the equations that define our final 2D magnetohydro- 
static equilibrium. The results from the numerical experiments 
are analysed closely in Sec.|4| whilst in Sec. |5]we comment on 
the nature of an analytical description for the equilibrium field. 
Finally, we conclude with a general overview of the problem in 
Sec.|6] 

2. Numerical scheme and initial setup 

For the numerical experiments studied in this paper, we have 
used Lare2D, a staggered Lagrangian-remap code with user con- 
trolled viscosity, that solves the full MHD equations, with the 
resistivity set to zero (see Arberetal. 2001). In order to cre- 
ate the initial magnetic field, a current-free hyperbolic X-point, 
A- - (x^ - y^)/2, is perturbed by squashing it in the vertical 
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y-direction by a given amount (1 - h) times the height of the 
original system, without introducing any plasma flow, such that 
the flux function of the initial state is given by 



(1) 



The squashing creates a uniform non-zero current density whose 
z-component is 

Mx,y,0)^^-1. (2) 

The squashing is characterised by the height of the box h, nor- 
malised to the original height, so that the dimensions of the ex- 
periment are L x hL. The initial plasma pressure po, is set to a 
constant everywhere, such that the initial system is not in equi- 
librium. 

The numerical domain is a 2D box with a uniform grid of 
1024x2048 points. The higher resolution in the y-direction is 
chosen to permit any current layer that may form to be as thin as 
possible, but still resolvable as far as possible across its width. 

Magnetic field lines are line-tied at the four boundaries and 
all components of the velocity are set to zero on the boundaries. 
The other quantities have their derivatives perpendicular to each 
of the boundaries set to zero. Hence, quantities that are con- 
served over the whole domain are total energy and total mass. 
Since the process is ideal (there is no diffusion to within the nu- 
merical limits), the field is frozen to the plasma, and mass in a 
single flux tube (or along a field line) must be conserved. The 
integrated current density is also conserved throughout the dy- 
namical evolution. This can be easily shown, taking the symme- 
try properties of the system into account. 



dt Js' lu dt Js 



V X B ds 



= ii<fiB.d, 

= -- (bB,dl = 0. 



(3) 



Here, S represents the surface of our numerical domain, with ds 
being a vector perpendicular to it, C is the the curve that encloses 
S , and dl is a vector tangent to C at each point. The integral of 
the tangential magnetic field, B,, along the curve C, which is the 
boundary of our box, equals zero for each pair of boundaries 
(left-right and top-bottom), because of symmetry. 

The numerical code uses the normalised MHD equations, 
where the normalised magnetic field, density and lengths, 

X = Lx , y - Ly , R - B„R , p = p„p , 

imply that the normaUsing constants for pressure, internal en- 
ergy and current density are 

Bl Bl . B„ 

Pn = — , e« = and ;„ = — . 

The subscripts n indicate the normalising constants, and the hat 
quantities are the dimensionless variables with which the code 
works. The expression for the plasma beta can be obtained from 
this normalization as 



B2 



In this paper, we will work with normalised quantities, but the 
hat is removed from the equations for simplicity. 



The (normalised) equations governing our MHD dynamical 
processes are 

^ + V ■ (pv) = , (4) 
ot 

d\ 

+ P(v ■ V)v = -Vd -h (V X B) X B -I- , (5) 

ot 



dp 

— -H V ■ = -yp^ - y + Hy , 
ot 



dB 

dt 



V X (v X B) 



(6) 
(7) 



where Fy and Hy are the terms for viscous force and heating, and 
internal energy is given by the ideal gas law, p = pe(y - 1), with 

r = 5/3. 

We have run a number of experiments with various heights 
ranging from h - 0.9 to h - 0.6, with the subsequent initial cur- 
rent densities ranging from jo - 0.23 to jo = 1.78, and various 
initial plasma pressures ranging from po = 1.0 to po - 0.125. 
Following a systematic study, we find that the function !F(A ) in 
equations ( fTOl i and ( fTTT i differs from one experiment to another, 
and is therefore determined by the initial plasma pressure and 
current density of the system. 

3. Magnetohydrostatic equilibrium in 2D 

In order to understand, and be able to give a description of, the fi- 
nal equilibrium state of our numerical experiments, we will now 
describe the basic equations of a two-dimensional MHS equilib- 
rium. 

In the absence of gravity, flows or frictional forces, the mo- 
mentum equation reduces to the following magnetohydrostatic 
(MHS) equation. 



jxB-Vp = . 



(8) 



where j and B are the electric current density vector and mag- 
netic field vector, respectively, and p is the plasma pressure. For 
our 2D problem, all the quantities depend on the horizontal and 
vertical coordinates, namely x and y. Making use of Ampere's 
law, V X B = /ioj, the 2D equihbrium can be characterised by the 
Grad-Shafranov equation. 



dp 

dA, 



1 



(9) 



which indicates that, in equilibrium, the plasma pressure is a 
function of the flux function, A^, where B = V X (0, 0, A,(x,3')), 
so such a state is uniquely determined by a function of the form 



(10) 



where f is an unknown function that is dependent on the initial 
conditions and evolution. Since, in 2D, the magnetic field lines 
are defined by contours of A,, the Grad-Shafranov equation tells 
us that in a MHS equihbrium, the plasma pressure is constant 
along field lines. 

In addition, it can easily be shown that V^A, = -j^, where 
j- is the z-component of the electric current density. Hence, j, 
should be also characterised by a unique function of A,, 



(11) 



when the system is in equilibrium, so that current density is also 
constant along field lines. 
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Fig. 1. Time evolution of the energies of the system, integrated 
over the whole two-dimensional box, for an experiment with 
jo = 1 -04 (h - 0.7) and po = 0.250. The magnetic, internal and 
total energies have been shifted on the y-axis by subtracting the 
values 0.238, 0.742 and 0.984 respectively, but their amplitudes 
are not to scale. 




Fig. 2. Two-dimensional contour plot of plasma pressure for the 
final equilibrium state for the sample experiment with h - 0.7 
and po = 0.250. White solid lines are the magnetic field lines as 
contours of the flux function A . 




Fig. 3. Two-dimensional contour plot (right) and surface (left) of electric current density for the final equilibrium state for the sample 
experiment with h = 0.7 and po - 0.250. White solid lines are the magnetic field lines as contours of the flux function A^. 



4. Numerical Results 

4.1. Energetics 

Initially, the only force acting on the system is the magnetic 
Lorentz force. This acts on the plasma by rising its kinetic en- 
ergy, which is zero initially. The velocities are then damped out 
through the viscous forces, and the physical MHD evolution al- 
lows the transfer of the kinetic energy into internal energy (heat), 
via the viscous heating term. This is a key process in the time 
evolution of the system, as it ensures the conservation of total 
energy, and hence, it will play a role in the final state. 

Fig. [T] shows the time evolution of the four energies of the 
system: kinetic, magnetic, internal and total, integrated over the 
whole two-dimensional box, as a function of the computation 
time normalised to the fast magnetoacoustic time, f/, defined as 
the time for a fast magnetoacoustic wave to reach the center of 
the box from the bottom or top boundary. 

The overall period of oscillation of the magnetic energy is 
about two fast magnetoacoustic time units, while the kinetic en- 



ergy, which is proportional to the velocity squared, has a pe- 
riod of oscillation of one fast magnetoacoustic time unit. The 
frequency ratio of these two oscillations is consistent with the 
Fourier analysis of the system, which involves a factor of two 
for the frequency associated to the internal energy. 

At the final equilibrium, where kinetic energy goes to zero, 
there has been an exchange from magnetic to internal energy. 
Magnetic energy is converted to kinetic energy through the mo- 
mentum equation, and then this kinetic energy is converted into 
internal energy (heating) through the viscous heating term. Total 
energy is conserved throughout the whole process. 

4.2. Final state 

In Figures |2] and |3] we show two-dimensional contour plots of 
the plasma pressure and electric current density, respectively, in 
the final state (defined here at f = 100, or equivalently, about 
290 fast-magnetoacoustic-wave transits), for an experiment with 
jo - 1 .04 (h = 0.7) and po - 0.250, together wit a three dimen- 
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(a)p„ = 0.750 (b)p„ = 0.625 (a)p„ = 0.750 <b)p„^ 0.625 




Fig. 4. Plots of electric current density across the width of the 
central current layer, for six different experiments, with h = 0.7, 
but different initial plasma pressures. 



Fig. 5. Plots of electric current density along the length of the 
main current sheet, for six different experiments, with h - 0.7, 
but with different initial plasma pressures. 




Fig. 6. The six plots in Fig. |5] are overplotted for comparison. The dimensions of the Green's potential solution are given by the 
dotted lines. Figure (b) is a zoom of (a) over a smaller range of current densities. In (b), the initial current density is overplotted 
(dashed). 



sional surface of the current density. Departing from an initial 
state containing an X-point with uniform pressure and current 
density, we get an equilibrium where the X-point has produced 
a thick current layer from which arms of enhanced current ex- 
tend along the curved separatrices (Fig.[3]l. The separatrices form 
cusp shapes at the two ends of the current layer. The plasma 
pressure is enhanced within the cusps (to the left and right of 
the current layer), and decreased in the regions outside the cusp, 
above and below the current layer (Fig.|2]i. Plasma pressure ap- 



pears to be constant along field lines. The electric current density 
is clearly not constant along the separatrices, but the surface cur- 
rent -current integrated over two given field lines at each side of 
the separatrix- is checked to be constant along them. 

We now look closely at the dimensions of the current layer 
for six experiments with jo = 1 -04 (h = 0.7) but different initial 
plasma pressure. In Fig. H] we show vertical cuts of the current 
density across the central current layer. The width of the cen- 
tral layer decreases for smaller plasma pressures, but remains 
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finite. In Fig. |5] we show horizontal cuts of the current density 
along the central current layer, for the same cases as in Fig.|4] As 
the pressure is decreased, the length of the central current layer 
extends further, and the current density becomes more concen- 
trated, developing a higher peak. This behaviour can be also ob- 
served if the initial plasma pressure is held fixed, and the height 
of the box is systematically decreased (i.e. the squashing is in- 
creased). Decreasing the initial plasma pressure has a similar ef- 
fect as increasing the initial current density, as the action of both 
is to make the Lorentz force dominate over the pressure force. 

When the initial plasma pressure is small (e.g. in Fig. |5]F), 
the current layer has a length that is many times longer that its 
width. We consider whether the current layer is approaching the 
form found in Green's current sheet solution. This is checked 
by comparing these plots with the le ngth of Green's sh eet (Fig. 
|6ll, given, by comparison with Bung ev & PriestI (Il995h . as a = 
0.25 -\fjo/7T, which for the sample experiment in Fig.|6] with jo - 
1.04 (h - 0.7), gives the value a ^ 0.144. 

In Fig. |6] the six horizontal cuts of Fig. |5] are overplotted, 
and the dimensions of Green's potential solution are marked. All 
the curves cross at the same points on x and y, namely {+a, jo), 
corresponding to the initial value of the current density in y, and 
the two ends of Green's current sheet in x. The main conclusion 
that may be extracted from these plots is that the field is in all 
cases very far from the potential solution, although the fact that 
all curves cross at the ends of Green's potential sheet seems to 
imply that Green's solution might be achieved (as far as we can 
get with the resolution) in the limit po — > 0. 

4.3. Current singularity 

The system has now reached a quasi-static equilibrium, and it is 
in an asymptotic regime where the current at the location of the 
null and along the separatrices slowly grows as time elapses. To 
evaluate the formation of a singularity we must look closer at 
the evolution of the current at the location of the null, namely, 
the "peak current". In Fig.|7]we show the time evolution of the 
peak current, j,„a^, for six different experiments with different 
initial plasma pressures. The first part of the evolution (from 
f = to f ~ 25) consists of the propagation of magnetoacoustic 
waves, bouncing back and forward between the boundaries and 
in the magnetic separatrices, which are quickly damped out by 
the viscous forces. After this relaxation has finished, the field is 
in equilibrium everywhere outside the separatrices, and the peak 
current at the null follows a very slow time evolution of the form 

= A(B + Ctf , (12) 

with D > 0. Both the maximum current, and the exponent D 
decrease as the plasma pressure is increased. That is, if a singu- 
larity is being formed, a large pl asma pressure hinders the for- 
mation of such (as discussed by ICraig & Litvinenkoll2005h . In 
any case, the formation of a singularity is a slow process that oc- 
curs in an infinite time, and, hence, it is impossible to reach, no 
matter how good the numerical resolution is. The key point here 
is that, once that the asymptotic regime has been achieved, any 
further grid improvement will only permit a longer run, but will 
not change the over all resu lts and conclusions. These results are 
in agreement with iKIappej (Il997h . who rigorously showed that 
for 2D ideal incompressible plasmas, a singularity cannot form 
in a finite time. In Fig [8] we show the evolution of the peak cur- 
rent for the same experiment but with different grid resolutions. 
It is shown that for both the value of the peak current, j^ax, and 
the growth rate, D, the effect of varying the grid resolution at this 
level is small. 



It is worth no ting that MHD re sistive instabilities such as 
the tearing mode ("Furt h et aDll973h may occur, halting the re- 
laxation and dissipating the current layer. In the idealised cases 
presented in the literature, the onset of this instability is inde- 
pendent of the resistivity and will occur, theoretically, for all val- 
ues of ?7, when the length-to-width ratio of the current layer has 
gone below 2n. It is found that our final quasi-equilibrium state 
is close to the condition given by this very simplified analyti- 
cal model. However, we have carefully checked that for all the 
results shown, no significant numerical reconnection has taken 
place at the location of the null, since the changes of the flux 
function at the origin remain below a relative value of 10"^. 

The study of such singularities is a complex topic, as they 
are, by definition, impossible to reach. A good way of ap- 
proaching this problem is by using mesh-refinement codes 
(e.g. Friedel et al. 199^ iGraue r & MarHani 1998) in which 
the grid resolution can be finely adjusted at certain locations. 
Nevertheless, these methods are not required for our purposes 
in this paper, as our high fixed resolution permits us to fully re- 
solve the current accumulations both at the null point and along 
(and across) the four separatrices (see Sec. 14. 4b . Moreover, it has 
been shown that the evolution of the singularity does not change 
its behaviour when different grid resolutions are used. 
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Fig. 7. Magnitude of the electric current density at the location 
of the null, as a function of time, for six different experiments, 
with h = 0.7, but with different initial plasma pressures. 
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Fig. 8. Magnitude of the electric current density at the location 
of the null, as a function of time, for the same experiment with 
h = 0.7 and po = 0.250, but with different grid resolutions. 
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(a) Horizontal cut 



(b) Vertical cut 
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Fig. 9. Electric current density across the width and along the length of the central current layer, for the same experiment with 
h = 0.7 and po = 0.375, at two different times, f - 100 (black diamonds) and t - 50 (dark red stars). These examples have the 
resolution of 1024 x 2048. 




Fig. 10. Electric current density across the width and along the length of the central current layer, for the same experiment with 
h - 0.7 and po = 0.375, at the time f = 50. Results from the high resolution run, 1024 x 2048 (black diamonds) are compared with 
the results of a half -resolution run, 512 x 1024 (dark red stars). 




Fig. 11. Residual forces along the central current layer, for the same experiment with h - 0.7 and po - 0.250, at (a) two different 
times, f = 100 (black diamonds) and f = 50 (dark red stars), and (b) with two different resolutions, 1024 x 2048 (black diamonds) 
and 512 X 1024 (dark red stars). 
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Fig. 12. Vertical cut of the current density at x = 0.1 for the 
experiment with h - 0.7 and po = 0.375, at two different times, 
f = 100 (black diamonds) and f = 50 (dark red stars), with the 
resolution of 1024 x 2048. The horizontal dotted lines at the top 
indicate the maximum currents in this cut. 
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Fig. 13. Test of Grad-Shafranov condition for the experiment 
with h = 0.7 and po = 0.375. Current density (black, y-axis 
on the left) and plasma pressure (blue, y-axis on the right) are 
plotted against the flux function A,, for every single point in the 
numerical domain. Positive values of Aj refer to inside of the 
cusp, while negative are outside the cusp. The functional form 
jzi^z) - jmaxi^ + mA,y" is ovcrplottcd here. 



4.4. Current layer 

The discussed on-going formation of an infinite time singularity 
suggests that small residual forces are acting at the location of 
the null, and hence, that the current layer is not in a perfect static 
equilibrium. In Fig. |9]we look at the dimensions of the central 
current layer at different times. Our high resolution experiments 
allow us to resolve the central layer both in length and width. It 
is shown here that both the length and the width of the layer get 
smaller as time elapses. 

Also, we check the dimensions of the current layer at equal 
times but with different grid resolutions (Fig. [TOb . finding that 
while the peak current varies slightly with resolution, the di- 
mensions of the current layer (at half the maximum value) do 
not vary. 

Finally, we look at the residual forces at the location of the 
null. These are checked to be zero everywhere in the domain 
except for the region very closed to the null that we show here. 
Fig. [TT]shows the x component of the total force, i.e. j x B - Vp, 
along the length of the central current layer, at different times 
and with different resolutions. It is shown that forces become 
smaller with time, but also they concentrate in a smaller region, 
meaning that the current layer gets shorter with time. 

As for the resolution, the regions in which the forces act 
(which we identify with the length of the current layer) stays 
the same, but the amplitude of the forces increases with reso- 
lution. This is a sign that these forces are not due to numerical 
reconnection, but are due to a singularity being formed. 

In conclusion, the system seems to be forming a singular- 
ity, but the process is slow enough so that we can describe it as 
a quasi-static equilibrium at a given time, in which the current 
layer at the location of the null has a finite length, a finite width 
and a finite amplitude. As time gets very large, the current den- 
sity layer gets smaller and thinner, but the peak value increases. 
It can be shown (Fig.fTSli that the current along separatrices also 
increases with time, as the length of the central layer decreases. 
This may suggest that this current layer is feeding current into 
the separatrices. 

Our aim now is to look for an analytical description of the 
field about the null in this quasi-static equilibrium at a given 
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Fig. 14. Logarithmic plot of the exponent n for both regions in- 
side and outside the cusp, as a function of the initial pressure, for 
a set of six experiments with the same squashing h - 0.7, after 
the same time has elapsed. The exponents follow a power law 
with a very similar decay rate for the two regions. The slopes are 
-0.63 and -0.6 respectively for inside and outside the cusp. 



time, that would precede reconnection in a non-ideal environ- 
ment. 



5. Analytical description of the field 

5.1. Grad-Shafranov equation 

A direct check on the validity of our equilibrium may be done 
by testing the behaviour of the pressure p with respect to the flux 
function A-, and also the consequences of the Grad-Shafranov 
condition, (|9]), which states that the current density must be 
also a unique fun c tion o f A,. We follow the approach given by 
IVekstein & Priest! (11993 ^ in an attempt to give a mathematical 
description of the field about the null in our quasi-static equilib- 
rium. They suggested the form jz{A^) = niA^", with n > 0, which 
is singular at the location of the null, where A, = 0. The problem 
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is that the state that we want to describe is not exactly singular 
at the null (even if it is slowly converging towards a singularity). 
Therefore, we try a function of the form 



jz(Az) = j„,ax(l+mA,y 



(13) 



where y„,a j is the value of the peak current at the null at a given 
time, and n > 0. This function is such that when mA, » 1, the 
field behaves as if it were singular at the origin, with the form 



(14) 



but when niA, «: 1, then j, = j„,ax- The value of m needs to 
be very large to counteract the small values of A, about the null. 
Note, that the separatrices do not satisfy the form given in Eq. 
[T3l as A; = but j, + j,nax along them. 

In Fig. [13] we represent every single point of the two- 
dimensional domain, for the plasma pressure and the current 
density against the flux function. It appears clear from this graph 
that the pressure is a unique function of A,, and so is most of the 
current density distribution. The biggest dispersion occurs when 
we approach A, - 0, which is the value on the separatrices, and 
at the X-point, where a singularity may be being formed, which, 
from the Grad-Shafranov equation, implies an infinite derivative 
of the pressure with respect to A,. 

On Fig.[T3] we have overplotted the best fit to Eq. ( fT3] ) to the 
current density, using i„ax - 8.5, and overplotted its derivative 
with respect to A- to compare to the plasma pressure. The best 
fit values, independently determined for the regions inside and 
outside the cusp, are: 

m; = 55 X 10"* 
n,- = 0.253 



m„ 
nn 



17 X lO"* 
0.294 , 



where the subscripts / and o denote inside and outside the cusp, 
respectively. It is clear that the form of the equilibrium is differ- 
ent for the regions inside and outside the cusp. These fits can be 
done for all the experiments with different initial plasma pres- 
sures or initial current densities. In Fig. [141 the exponents n are 
plotted for different initial plasma pressures for both inside and 
outside the cusps. These exponents follow in both cases a nega- 
tive power law with a very similar slope, suggesting, once more, 
that an increase in plasma pressure acts against the formation of 
a singular current density. 

5.2. Poloidal flux function 



IVainshtein (1990) proposed a description of the field about spe- 
cial magnetic points in two-dimensions, such as cusp points at 
the ends of thin current sheets, seeking a solution at small r of 
the form 



A,{r,e) = r«>gi(0) + r'^'-giie) + ... + r^'gM 



(15) 



with 1 < ai < q;2 < ■■■ < Cn, thus avoiding A; —> oo when r — > 0. 
Here, r is the radial coordinate whose origin is at the location of 
the null, and 6 is the angular coordinate which is zero at the main 
axis of the current layer (i.e. jc-axis). 

We suggest a poloidal flux function of the form 



A,(r,6) ^Ao{r)+Ai(r,6) , 



(16) 
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Fig. 15. Analytical prediction of the flux function, compared to 
numerical equihbrium state, for the experiment with h - 0.7 and 
Pi) — 0.375. In (a) we show the numerical magnetic field lines 
(solid black) over a contour map of the current density, with the 
analytical contours of A. for inside the cusp (red dashed) and 
outside the cusp (blue dashed). In (b) and (c) we show the rel- 
ative error of the analytical prediction from the numerical state 
for a horizontal cut (inside cusp) and a vertical cut (outside cusp) 
respectively. 



where 

Ao(r) 
Ai(r,e) 



ar'' , 

br'' cos 26 . 



(17) 
(18) 



This form can be compared to Eq. ( [T3b . using V^Aj. - -j,. From 
Equations (ITSI l. ([TTl l and (ITSl l. we get 



V^A, 



ap^r'''^ -H biq- 



cos 261 



(19) 



Now, assuming the singular form for the flux function given by 
Eq. ([T4l l. combined with Eq. (ITSI l. we obtain three relations be- 
tween the parameters a, p, q and the parameters m and n, as 
follows. 



I + n 



m"p-^ 



l/in+l) 



(20) 
(21) 
(22) 



Therefore, the only free parameters of our analytical fits are 
(m, n, b). In Fig. [TSh . we show the contours of the numerical 
Aj, compared to the contours of our analytical form of A.(r, 9). 
The fit is fairly good save at the origin and near the separatrices. 
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where the assumptions above are not valid, since the field is not 
singular there. Figures [T5b and [TSb show the relative deviation 
of the analytical solution from the numerical in horizontal and 
vertical cuts at the center of the box respectively. The fits along 
these cuts are very good, with a very small error of the order of 
10"^ - 10"'*. The parameters of the fit are summarised in Table 
[U for the experiment with h = 0.7 and po - 0.375, and with 
jmax = 8.5. The regions inside and outside the cusp are denoted 
with / and o respectively. 

These parameters depend on the initial plasma pressure and 
the initial current density, but also on the time we choose for 
our quasi-static equilibrium (i.e. they depend on the choice of 
the peak current). However, they provide a qualitative result for 
an non-force-free X-point equilibrium, shortly before magnetic 
reconnection occurs. 



Table 1. Parameters of the analytical solutions for the regions 
inside (/) and outside (o) the cusp. 





111 


n 


h 


a 


P 


q 




55 X 10-' 


0.253 


0.60 


-0.18 


1.60 


1.83 


o 


17 X 10-* 


0.294 


0.66 


-0.17 


1.55 


1.82 



6. Summary and conclusions 

We have presented evidence that the ideal MHD relaxation of 
two-dimensional magnetic X-points embedded in non-zero beta 
plasmas leads to a non-force-free quasi-static state that depends 
highly on the initial values of the background plasma pressure 
and electric current density. In this state, the system is in equi- 
librium everywhere save at the null and along the separatrices. 
Hence, the evolution of the system can be divided in two parts: 
(1) the quick viscous relaxation of the whole system, plus (2) 
an asymptotic regime with a very slow evolution of the current 
density at the null and along the four separatrices. 

The slow increase of the current density at the location of 
the null is associated with the evolution to a singularity in an 
infinite amount of time. This slow evolution at the null makes the 
current accumulation region shorter and thinner as time elapses, 
concentrating the current at the exact location of the null, and 
also along the four separatrices. It should be noted that the large 
numerical resolution used here permits us to resolve the current 
layer both in width and length. However, note also that a singular 
value of the current will always be impossible to reach because 
(i) the singularity does not form in a finite time and (ii) any grid, 
no matter how fine, will eventually result in a numerical diffusion 
of the current. 

Our overall re sults are completely_ different to those 

of [R astatter et al.^ 



991 . ICraig & Litv inenkd (l2005h and 



IPontin & Crai^ (I2005h for two main reasons: First, as time gets 
large, the system is not trying to collapse to a current sheet with 
a finite length, but it is concentrating all the current at the loca- 
tion of the null and along the separatrices. Second, this infinite 
current state is not achievable, and therefore, we must describe a 
final state in a quasi-static equilibrium at a given time, in which 
the current density has accumulated at the location of the null 
in a layer with finite length and width, that we are able to re- 
solve in our numerical experiments. These results suggest that 
the idea of an infinitesimally thin current sheet must be then put 
aside when describing the MHD evolution of the system. The fi- 
nal state described in this paper may be identified with a moment 



shortly before magnetic reconnection takes place in a consistent 
non-force-free two-dimensional model. 

We propose a functional form of the quasi-static equilibrium 
izi^z) - 7max(l + niA^)"" , and we derive a form for the flux 
function of A.(r, 6) - ar^ + br^ cos 20, where the parameters 
a, p and q are well defined functions of m and n. It is found 
that all these parameters are different from the regions inside and 
outside the cusp, and hence, the final equilibrium is different, 
and the system approaches the singularity in a different way for 
positive and negative values of A^, i.e. inside and outside the 
cusp, respectively. Also, the final equilibrium depends strongly 
on the initial plasma pressure and current density. None of the 
above expressions are valid along the four separatrices, where 
A, = but 7j jmax- Hence, extra considerations may be needed 
for a full description in those regions. 

It is found, both numerically and analytically, that an in- 
crease in the initial plasma pressure, and hence, an increase in 
the plasma beta, acts on the field by reducing both the strength 
of t he singularity and its grow th rate with time. This suggests 
(see 'Craig & Litvinenkdl2005h that plasma pressure acts against 
fast reconnection. 

Even if this paper describes only a qualitative analysis of 
a 2D hydromagnetic equilibrium, it describes a fair approxi- 
mation of the behaviour of the final state as the values of the 
initial plasma pressure and current density are varied. These 
two-dimensional contexts are of high relevance for systems with 
translational or rotational symmetries, and their study is useful 
for some astrophysical environments which can be well approx- 
imated by these properties of symmetry. 

In a follow-up paper, we will evaluate current accumulations 
in three-dimensional equilibria which contain 3D magnetic null 
points. The characteristics of these environments are going to 
be completely different to the two-dimensional case, and the 
dynamical evolutions are less restrictive in the sense that the 
plasma has freedom to move in all three spatial directions. 
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